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Using transfer matrix and finite-size scaling methods, we study the thermodynamic behavior of a 
lattice gas with two kinds of particles on the square lattice. Only excluded volume interactions are 
considered, so that the model is athermal. Large particles exclude the site they occupy and its four 
hrst neighbors, while small particles exclude only their site. Two thermodynamic phases are found: 
a disordered phase where large particles occupy both sublattices with the same probability and an 
ordered phase where one of the two sublattices is preferentially occupied by them. The transition 
between these phases is continuous at small concentrations of the small particles and discontinuous at 
larger concentrations, both transitions are separated by a tricritical point. Estimates of the central 
charge suggest that the critical line is in the Ising universality class, while the tricritical point 
has tricritical Ising (Blume-Emery-Griffiths) exponents. The isobaric curves of the total density as 
functions of the fugacity of small or large particles display a minimum in the disordered phase. 
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I. INTRODUCTION 

Models for fluids with only repulsive interactions have 
been studied in the literature for quite a long time [l|, 
both on a lattice and in the continuum. Although at¬ 
tractive interactions are essential to produce liquid-gas 
transitions, models with only repulsive interactions may 
show transitions which resemble the melting of a solid 
phase. If the repulsive interactions are of the hard-core 
excluded volume type, all allowed configurations of the 
system have the same energy and thus the model is ather¬ 
mal. The continuous versions of these models are known 
as hard sphere models, and a fluid-solid transition was 
found in their phase diagram. It is worth recalling that 
in the seminal work by Metropolis et ah, where Monte 
Carlo simulations were introduced, the hard sphere gas 
was studied [I|, 0 • 

Athermal lattice models with excluded volume interac¬ 
tions have also been widely studied before. In these mod¬ 
els the localization of particles is constrained to sites of 
a lattice. A model where a particle excludes others from 
only its site corresponds to an Ising lattice gas without 
the interaction term, and no phase transition is found. If 
a particle placed on a site of the square lattice excludes 
other particles from its four first neighbor sites, at low 
densities the sites of both sublattices are equally occu¬ 
pied, but as the density is increased one of the sublattices 
will be preferentially occupied by the particles. One may 
associate the disordered low density phase to a fluid and 
the ordered phase to a solid, so that the transition may 
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be seen as a melting of a solid phase. This model has 
been thoroughly investigated by a variety of techniques 
0. On the square lattice, the universality class of the 
transition was object of some discussion in the literature 
0. Accurate estimates of the thermodynamic behavior 
of the model using transfer matrix and finite-size scal¬ 
ing techniques show a continuous transition between the 
disordered and the ordered phase, which is in the Ising 
universality class as one should expect from sym¬ 
metry considerations. This model may be generalized 
by increasing the range of the excluded volume interac¬ 
tions, and discontinuous phase transitions are found if 
this range is large enough. We refer to recent simula- 
tional investigations of this family of models where com¬ 
prehensive surveys of the literature may be found in. 
Although the transfer matrix formalism may lead to pre¬ 
cise estimates for the critical behavior of such athermal 
models, Monte Carlo methods with cluster algorithms 
may also furnish good results. As examples, we mention 
the IsiM lattice gas with first neighbor exclusion on the 
cubic U and triangular @ lattices. 

Here we consider another generalization of the lat¬ 
tice gas with first neighbor exclusion, introducing also 
small particles which exclude only the site they occupy. 
This mixed lattice gas was studied, using series expan¬ 
sion techniques on the square lattice, by Poland 0- He 
found evidences for a tricritical point in the phase dia¬ 
gram of the model: while for low densities of small par¬ 
ticles the transition is continuous, it becomes discontin¬ 
uous as the density is increased. A slight modification of 
this model, where the large particles occupy elementary 
squares of the lattice and the small particles are located 
on the center of edges, has been exactly solved in the 
grand canonical formalism. In a particular case of this 
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model, when the fugacities of the small (zi) and large 
(Z 2 ) particles obey the relation Z 2 = (1 + zi)^, Frenkel 
and Louis were able to show that it may be mapped on 
the Ising model with vanishing magnetic field and thus 
its solution in two dimensions is known 0. It should be 
mentioned that the fact that a model very similar to the 
original can be mapped on the Ising model at a partic¬ 
ular point of the critical line is a strong additional evi¬ 
dence that the transition is in the Ising universality class. 
In particular, for Zi = 0 the Frenkel-Louis model corre¬ 
sponds to the model with only large particles when their 
fugacity is unitary. A variant of this model was proposed 
and studied by Lin and Taylor [I^. In this model, the 
small particles are triangles such that one of their sides 
occupies an edge of the lattice and the opposite vertex is 
located at the center of an elementary square, so that up 
to four small particles may be placed in an elementary 
square of the lattice. If finite interaction energies between 
the large (square) particles on first-neighbor sites and be¬ 
tween large (square) and small (triangle) particles which 
share a lattice edge are introduced, this model shows a 
lower critical solution point. We also notice that a deco¬ 
rated lattice model for hydrogen-bonded mixtures which 
can be mapped on the three-dimensional Ising model and 
which shows both upper and lower critical solution points 
has been studied by Wheeler and Anderson (l^ . 

More recently, the solution in the grand-canonical for¬ 
malism of the model on a Bethe lattice with arbitrary 
coordination number was obtained [l3| . The general fea¬ 
tures of the phase diagram are consistent with the find¬ 
ings of Poland 0- In the parameter space defined by 
the two activities Z 2 and zi of large and small particles, 
respectively, a continuous transition is found for small 
values of zi, and it becomes discontinuous if zi is suffi¬ 
ciently large. Since the slope of the critical line is neg¬ 
ative at low values of zi, becoming positive at higher 
values, a re-entrant behavior of the transition is seen in 
this region of the phase diagram. At even higher values 
of the activity of small particles, the transition becomes 
discontinuous and thus a tricritical point is found. An¬ 
other interesting feature of the model is that the isobaric 
curves, where the total density of particles is considered 
as a function of the density or the activity of small par¬ 
ticles, shows a minimum in the fluid phase. We recall 
that perhaps the most studied situation in nature where 
a density anomaly is found happens in water close to the 
freezing point, where a maximum of the isobaric curves 
of the density as a function of the temperature is seen. 
The Bethe lattice solution of the mixed lattice gas model 
leads to a minimum of the density as a function of an¬ 
other held-like thermodynamic variable: the activity of 
small particles. Although the physical situation in wa¬ 
ter is of course quite distinct from the model we study 
here, it is noteworthy that several recent studies of sim¬ 
ple effective models for water suggest that the density 
anomaly may be due to effective interparticle interaction 
potentials with two length scales [T^, a feature which is 
also present in the mixed lattice gas model. 
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FIG. 1. (color on line) Strip of width L = 4, with periodic 
(dashed-red) and helical (dotted-blue) boundary conditions. 
The lower and upper rows of sites determine the states of 
the transfer matrix. They are (0, 2, 0,1) and (0, 0,0, 1), re- 
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spectively. The corresponding element is Z 2 ' zi. The par¬ 
tially overlapping dot-dashed boxes are the states ((2, 0,1, 0) 
and (0,1,0, 0))for helical boundary conditions, and the cor¬ 
responding matrix element is 1, since a single empty site is 
added. 


Through transfer matrix and finite-size scaling calcula¬ 
tions, we show here that the thermodynamic behavior of 
the model on the square lattice is qualitatively the same 
found in the Bethe lattice solution. Our results suggest 
that the critical line is in Ising universality class, simi¬ 
larly to the model with large particles only, and that the 
tricritical point belongs to the tricritical Ising (Blume- 
Emery-Griffiths - BEG) class. 

In Sec. im we define the model more precisely and 
describe its transfer matrix solution on strips of finite 
widths, for both periodic and helical boundary condi¬ 
tions. We also discuss the factorization of this transfer 
matrix in a product of sparse matrices, which reduces the 
computational effort to handle them. Our results for the 
phase transitions and other thermodynamic properties of 
the model in the two-dimensional limit may be found in 
Sec. mi Section m is devoted to final discussions and 
conclusion. 


II. DEFINITION OF THE MIXED LATTICE 
GAS MODEL AND ITS TRANSFER MATRIX 

We will use the transfer matrix formalism to study 
the model with both particles placed on a square lattice. 
This will be accomplished by solving the model on strips 
of finite width L. For periodic boundary conditions, the 
states of the transfer matrix will be defined by the config¬ 
uration of the L lattice sites in the same row of a cylinder. 
As an example, for L = 4 we have a total of 26 states, but 
this number reduces to 9 if the rotation symmetry is con¬ 
sidered. These states, with their multiplicity indicated 
between curly brackets, are: (0, 0,0, 0) — {!}, (0, 0,0,1) — 
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{4}, (0, 0,0, 2) - {4}, (0,0,1,1) - {4}, (0,1, 0,1) - {2}, 
(0,1,0, 2) - {4}, (0,1,1,1) - {4}, (0, 2,0, 2) - {2}, and 
(1,1,1,!) — {!}. In our notation, 0 represents an empty 
site, while 1 and 2 correspond to sites occupied by small 
and large particles, respectively. To avoid frustration in 
the solid phase, we restrict ourselves to even widths. The 
transfer matrix is obtained considering two adjacent rows 
of L sites in a particular configuration of particles and 
checking if the excluded volume interactions are satis¬ 
fied. If this is the case, the corresponding element of the 
transfer matrix will be given by Tj,(i, j) = ^ 

where ni and ni are the total numbers of small and large 
particles in both rows, respectively, and K is the number 
of ways (considering the excluded volume) of placing the 
row 2 over the row i, by keeping % fixed and rotating j. 
Defining the matrix elements in this way assures that the 
transfer matrix is hermitian. For L = 4, some elements of 
the transfer matrix, with the states ordered as presented 
above, are given by: 

Ti{l,2) = T4(2,2)=4zi, (1) 

r4(3,9) = 0, T4(5,8) = zi22- 

Besides adopting periodic boundary conditions, we 
also considered helical boundary conditions, where the 
rightmost site of a row is linked to the leftmost site of 
the row above. Both boundary conditions are illustrated 
in Fig. [I] As happens for periodic boundary conditions, 
for helical boundary conditions, the states are determined 
by the configuration of two sets of L sites, where the sec¬ 
ond one is shifted by one site only to the right of the first 
set, as is also shown in Fig. [TJ Thus, while for periodic 
boundary conditions at each application of the transfer 
matrix L new sites are added to the lattice, for in the he¬ 
lical case a single site is added. This has the advantage 
of leading to a sparser matrix, since L — 1 sites of both 
configurations are coincident. In Fig. [l] two successive 
site sets, used to define the transfer matrix elements for 
helical boundary conditions, are enclosed by dashed and 
dotted rectangles. 


A. Number of states 

The numbers of states for lattices of successive widths 
are related by simple recursion relations. Let us call 
the number of states in a strip of width L and 
free boundary conditions, such that the first site is in the 
configuration ai and the last site in the configuration ai,. 
Since (7^ = 0,1, 2, we have a total of 9 of these numbers of 
configurations, where 6 are independent due to reflection 
symmetry. The total number of configurations will be 
= So-i aL boundary conditions. For 

periodic boundary conditions, the occupancy of the sites 
1 and L should obey the excluded volume constraint, so 
that the number of states in this case will be given by 
= Scti aL ’ where the prime restricts the sum 

to terms such that cti + ctl < 2. The numbers of states 


obey a set of nine linear recursion relations: 

(2a) 
(2b) 
(2c) 

Iterating these recursion relations with the starting val¬ 
ues niWi = ^ 0 - 1 ,(TL j may obtain the numbers of states 
for increasing values of the widths for both boundary 
conditions, which are shown in Tab. ID The states for he¬ 
lical boundary conditions coincide with the ones for free 
boundary conditions. It is easy to reduce these recursion 
relations to a single one: 
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with the remaining numbers of states being given by 
arid The 

^01,0 — ”- 01,1 ^01,i> ana 2 — ^01,1 ’^oi,! ■ 

characteristic polynomial associated to the linear recur¬ 
sion relation (jS]) will be: — 2r^ — r + 1 = 0, which has 

three real roots. Although of course the roots may be 
found analytically, we will not give the expressions here 
since they are rather long. The approximate numeri¬ 
cal values are ri = 2.2469797, r 2 = —0.80193776, and 
ra = 0.55495811. Therefore, in general, we have that: 


"-o-i,! ~ ^01^1 


B, 


01'2 


+ a 


01 ' 3 ) 


(4) 


where the coefficients are determined by the initial con¬ 
ditions. To obtain the actual numbers of states for fi¬ 
nite widths, it is more practical to iterate the recur¬ 
sion relations directly, but explicit expressions such at 
o are useful to obtain the asymptotic number of states 
for L 1, which is dominated by the leading root of 
the characteristic equation, r; rf Ri 2.2469797^. 
The coefficient of this asymptotic behavior will be dif¬ 
ferent for free and periodic boundary conditions. This 
may be compared with the result for the model with 
large particles only, where we have, for periodic bound¬ 
ary conditions that the numbers of states are 
given by a Fibonacci sequence starting with 1 and 3 
ni{L) = Fi^ 3 {L) r ; [(1 -|- y/5)/2]^ r ; 1.6180-^. As ex¬ 
pected, the number of states increases much faster with 
the width when small particles are present. 

In order to reduce the amount of memory and com¬ 
puter time demanded, for periodic boundary conditions, 
we use a method proposed by Todo and Suzuki [l^ that 
consists in decomposing the matrix in L + 1 sparse 
matrices, so that Tl = Here, 

instead of adding a complete row in the strip, each ma- 
trix , with i = 2,...,L acts adding a new lattice 
site to the row. The matrix creates a site z = 1 
in a new row h, just above the site z = 2 from the row 
h — 1, accounting for the NN exclusion between particles 
in these sites. Then, ’ adds the site z = 2 in the row 
h (above the site [3, h — 1]), imposing the exclusion with 
particles in the site below it and the site z = 1 in the 
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L 

Uf 

Up 

1 

3 

2 

2 

6 

6 

3 

14 

11 

4 

31 

26 

5 

70 

57 

6 

157 

129 

7 

353 

289 

8 

793 

650 

9 

1782 

1460 

10 

4004 

3281 

11 

8997 

7372 

12 

20216 

16565 

13 

45425 

37221 

14 

102069 

83635 

15 

229347 

187926 

16 

515338 

422266 

17 

1157954 

948823 

18 

2601899 

2131986 

19 

5846414 

4790529 

20 

13136773 

10764221 


TABLE I. Numbers of states of the transfer matrix for strips 
of width L for free (n/) and periodic (rip) boundary condi¬ 
tions. 


row h. The other matrices up to work in the same 
way. Finally, the matrix imposes the NN exclu¬ 

sion between particles in the sites L and 1, and translate 
the labels i —>■ i — 1, completing the formation of the 
row h. More details may be found in Q, where this 
method was applied to the case of large particles only. 
As an example of the power of this method, the num¬ 
ber of non-vanishing terms in the symmetric matrix for 
L = 14 is Ms = 31895812, whereas the sum of such terms 
in the fifteen sparse matrices is Msp = 6706321, the ratio 
r = Mg/Mgp is r ft! 4.7. For L = 16 this ratio is r ft 15.6 
and for L = 18 it is r ft 56.8. This method allows us to 
work with strips of widths up to L = 18 in computers 
with 16GB of RAM. The numerical method we used to 
find the two dominant eigenvalues of the transfer matrix 
was a variant of the power method till 


After performing some numerical calculations for pe¬ 
riodic and helical boundary conditions, we realized that 
the first, combined with the factorization of the transfer 
matrix, are more efficient, allowing us to handle strips 
of larger widths with the same computational effort. We 
will, therefore, present results for periodic boundary con¬ 
ditions only. Notice that to perform the factorization of 
the transfer matrix, no block-diagonalization using the 
symmetries is done. 


III. RESULTS FOR PERIODIC BOUNDARY 
CONDITIONS 


A. Phenomenological renormalization 


The inverse correlation length of the model for a strip 
of width L, is given by: 


where Xi is the largest eigenvalue of the correspond¬ 
ing transfer matrix. From previous results [id . Il4| . we 
may expect to find continuous and discontinuous transi¬ 
tion lines between ordered and disordered phases. Phe¬ 
nomenological renormalization [T^ states that the esti¬ 
mates for these transitions can be obtained from the fixed 
point of the recursion relation: 


Cl _ Cl' 
L ~ U' 


( 6 ) 


As discussed above, due to the symmetry of the crys¬ 
talline phase, only even strip sizes must be considered, 
so that L' = L -I- 2 is used. The transition lines are 
shown in Fig. [2 where a very small dependence on the 
size L is observed. For a fixed value of zi, we may expect 
this dependence to be: 


• 22 (A) — Z2c + UlA -f 0,2^ . (7) 

where Z 2 c and the amplitudes ai should be functions of 
zi , while the exponents are expected to be constant along 
the critical line. For zi = 0, Guo and Blote @ found that 
Xi = 2 + i and, in fact, for small zi, we have verified that 
xi ft 3 and X 2 ft! 4. Then, assuming that these expo¬ 
nents hold for the whole transition curve, we obtain the 
values of Z 2 c from Eq. ©■ The resulting transition line 
(TL) is also shown in Fig. [51 but no clear difference is 
observed from the ones for a given L at this scale, due 
to the small finite-size corrections. Indeed, the maximal 
difference between the extrapolated TL and the line for 
the pair (16,18) is ft 0.1%, which gives us an idea of 
the error bars in these estimates. In accordance with the 
Bethe lattice solution of the model [l^, the critical line 
initially has a negative slope that approaches —1 when zi 
approaches zero. Increasing zi, this curve reaches a min¬ 
imum, located at zi ft 0.3805 and Z2 ft 3.5976 - larger 
than the one found in the Bethe lattice (zi ft 0.26 and 
Z 2 ft 1.55) - and then starts to increase. Finally, we ex¬ 
pect the curve to reach the tricritical point, beyond which 
the transition becomes discontinuous. The estimated tri¬ 
critical point is also shown in the figure, and the methods 
we used to obtain this estimate will be discussed below. 

The estimates for the location of the tricritical (TC) 
point using results for finite-size strips may be obtained 

from the condition ~ ^ ~ (T^ . The values 

L/ — 2 ij L/ -\- 2 

of and z'^'^ found for different L's are show in Tab. 
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FIG. 2. (Color on line) Finite-size and extrapolated transition 
lines, calculated from the fixed point of the phenomenological 
renormalization transformation, considering the pair of strip 
sizes L,L' = L -\-2. In the main plot, the differences between 
the estimates are almost not visible. In the inset the region 
close to the minimum is amplified, so that the finite-size ef¬ 
fects are seen. For a given value of z\, as the width of the 
strip increases the estimates for the critical value of Z 2 de¬ 
crease. The estimated tricritical point, discussed below, is 
represented by the black circle. 


L 

Zl 

Z2 

6 

2.879392 

11.838332 

8 

3.055064 

12.967494 

10 

3.194821 

13.918518 

12 

3.286902 

14.570770 

14 

3.348568 

15.019111 

16 

3.391375 

15.335858 

00 

3.559(2) 

16.60(4) 


TABLE II. Estimates for the location of the tricritical point, 

using the condition ^ — ^ = .§^+2. 

^ L-2 L L+2 


im Assuming that 

zf^{L) = zf^ + + ..., ( 8 ) 

where z = I, 2, we may find the exponents t/ip considering 
that bij = 0 for j > I and performing a three-point fit. 
Namely, we solve for zf^ , and in 

zf^{L + l) = zf'"+bi^i{L + l) ( 9 ) 

with Z = 0, — 2 and 2. The values of are not constant, 
but converge toward asymptotic values as L increases, as 
shown in Fig. [3l Due to finite-size effects, for L = 8, 
negative exponents are found, which do not make sense. 
Extrapolating the exponents for L > 8, we find the same 


FIG. 3. (Color on line) Exponents of finite-size corrections 
(Eq. ((8|) for the tricritical point activities estimated from 
phenomenological renormalization. In the main plot yip and 
y 2 ,i are extrapolated with A = 3.38 and A = 2.98, respec¬ 
tively. The behaviors of yip and y 2,2 are shown in the inset. 


asymptotic value ypi ~ 1.75 for both i = 1 and i = 2 
(see Fig. Therefore, we will assume that y^i = 7/4. 
Inserting this exponent in Eq. m and assuming that 
bij = 0 for j > 2 , the value of y^p can be obtained 
from a four-point fit (for the unknowns zf^, 6 pi, 6^2 
and yip)- These exponents are shown in the inset of Fig. 
|3l Unfortunately, they do not present a monotonic be¬ 
havior, possibly due to the small strip widths we could 
handle, so we cannot extrapolate them, but we may see 
that yip ^ 5 and y 2 p ~ 4. Considering different expo¬ 
nents in the ranges 5 ^ yip ^ 7 and 4 ^ y 2 p ^ 6 , we 
obtain several estimates of the tricritical point, leading 
to zf^ = 3.557(3) and Z 2 ^ = 16.60(2). These values and 
their respective error bars are also calculated considering 
different L’s in the fit (not only the three largest strip 
sizes). We notice that these values of the activities at the 
tricritical point are much larger than the ones found on 
the Bethe lattice with coordination q = A {zf^ = 1.16956 
and Z 2 ^ = 3.02938) [T^. It is, indeed, expected that 
mean-field approximations underestimate the critical ac¬ 
tivities. 

Beyond the tricritical point the transition becomes dis¬ 
continuous and, at first sight, we should not extrapolate 
the coexistence line in the same way we did for the crit¬ 
ical line (below the TC point). Actually, from finite-size 
scaling, we expect exponential corrections at coexistence 
[ 2 ^, instead of the power law behavior (e. g., Eq. ©) 
at criticality, namely 

^2{L) = Z2c + fllC -f . . . (10) 

However, performing a three point fit considering this 
correction form, we have found a coexistence curve whose 
maximal difference from the one obtained as above (using 
Eq. (O) is smaller than 0.02%. Thus, both approaches 
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From a three-point extrapolation considering these expo¬ 
nents, we obtain = 3.5541(3) and xj*" = 16.5827(3), 
which agree with the estimates from phenomenological 
renormalization shown above, but are more accurate. 
These values will hereafter be used as the location of 
the TC point. 

TableHnialso shows the central charge at the maximum 
as a function of L, where a non-monotonic behavior can 
be observed, preventing us to extrapolate this quantity 
in a reliable way. Nevertheless, the values in Tab. IIIII 
are very close to c = 7/10 and they are approaching this 
limit for the larger widths. 

C. Densities of particles 


FIG. 4. (Color on line) Central charge c as a function of 
zi, calculated from f{L') — f{L) along the finite-size critical 
line for the pairs L, L' indicated. The dashed horizontal lines 
indicate the Ising critical (bottom) and tricritical (top) values. 


The densities of particles for the fluid (F) and 
(S') phases can be obtained from 


Pf (L) = Zi 


^z^ 


and pf (T) = Zi 


df{L) 

dzi 


solid 

(13) 


lead to very similar coexistence lines in the thermody¬ 
namic limit. 

B. Conformal anomaly 

The free energy per site of the system for a strip of 
width L is given by 

/(L) = ilnAi. (11) 

Conformal invariance theory states that at criticality 


with 1 = 1,2, and f^{L) and f^{L) being the free ener¬ 
gies of the system in fluid and solid phases, respectively. 
We calculate these densities along the finite-size transi¬ 
tion lines (TL) in the following way: along the TL ob¬ 
tained for the pair (L — 2, L), we determine the densities 
for the strip width L. The resulting densities curves are 
displayed in Fig. [SK- For small zi (and pi - i. e., in 
the critical line) a very small finite-size dependence is 
observed in the curves for the F phase, while for large 
Pi (i. e, in the coexistence line) some L-dependence ap¬ 
pears. Assuming finite-size corrections in the form 

Pi{Li) = Pi -\- o-ipL -|- ai^2L -b ..., (14) 


f{L) ^ /oo + ^, (12) 

where c is the central charge. Figured shows this quan¬ 
tity as a function of Zi along the finite-size critical lines 
and we can see a crossover from c ~ 0.5 at small zi to 
c ~ 0.7 near the TC point. Moreover, as the widths of the 
pair of strips become larger, the estimates for c approach 
c = 1/2 at low values of Zi and the crossover is steeper. 
Therefore, in the thermodynamic limit we expect that 
c = 1/2 for zi < zf^ and c = 7/10 at Zi = in good 
agreement with the critical and tricritical Ising univer¬ 
sality classes, respectively, similar to what is found, for 
example, for the BEG model [2l|. 

The central charge curve, for a given pair (L, L'), has 
a maximum close to c = 7/10 that should be located ex¬ 
actly at the tricritical point when L —^ oo. Thus, we may 
dc 

use the condition —— = 0 as an alternative estimate of 

UZi 

this point. The values obtained in this way are shown 
in Tab. IIIII Performing extrapolations similarly to what 
was done in the previous subsection, we find pi^i ~ 2.37. 
Moreover, exponents yi ^2 ~ 4.5 and j/ 2,2 ~ 9.9 are found. 


from three and four point fits, we obtain the exponents 
Vi^i ~ 2.8 and Vi ^2 ~ 4. Then, using these exponents 
to extrapolate the transition line for the fluid phase, we 
find the extrapolated curves shown in Fig. [51 d. It is worth 
noticing that as zi and, consequently, pi goes to zero (on 


iL,L') 

Zl 

Z2 

Cmax 

(6,8) 

3.461580 

15.844423 

0.700552 

(8,10) 

3.503228 

16.172684 

0.702456 

(10,12) 

3.522518 

16.327104 

0.702901 

(12,14) 

3.532852 

16.410446 

0.702867 

(14,16) 

3.538962 

16.459910 

0.702677 

(16,18) 

3.542843 

16.491393 

0.702447 

00 

3.5541(3) 

16.5827(3) 



TABLE III. Tricritical point estimates from the maxima of 
the central charge curves (Fig. 2)) and the maximum values 
of the central charge Cmax ■ 
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L 

pi 

p2 

10 

0.3780(5) 

0.2038(3) 

12 

0.3780(7) 

0.2037(4) 

14 

0.3780(9) 

0.2037(5) 

16 

0.378(1) 

0.2037(7) 

18 

0.378(1) 

0.2037(8) 

00 

0.378(1) 

0.2037(8) 


TABLE IV. Densities of small pi and large p 2 particles cal¬ 
culated at the tricritical point. 


the critical line) the density of large particles approaches 
the value p 2 = 0.367742... estimated by Guo and Blote 
0 for pi = 0. 

Along the density curves for the solid phase, we may 
observe a re-entrant behavior, with the region below and 
above the turning point related to the critical and coex¬ 
istence lines, respectively. At odds with the fluid phase, 
in the coexistence region (large P 2 ) negligible finite-size 
effects are observed, while strong corrections appears for 
small p 2 (in the critical curve). In this last region, below 
the turning point, if we assume that corrections are given 
by Eq. (HI, we find the exponents ~ 1 and Vi ^2 ~ 3, 
which lead to extrapolated critical curves in very nice 
agreement with the ones for the fluid phase. This indi¬ 
cates that below the TC point we have in fact a critical 
situation where the densities of the two phases must be 
equal. Moreover, this confirms that the correction form 
assumed (Eq. (fTdll l with the exponents associated is in 
fact appropriate. 

Although the Eq. (fHl) still works on the coexistence 
line for the solid phase (large P 2 ), the exponents change 
to Vi^i ~ 0.3 and ~ 1.3. The extrapolated coexistence 
curves with these exponents are also shown in Fig. [SId. 

The densities of particles for strips of size L calculated 
at the extrapolated tricritical point are shown in Tab. El 
It is interesting that, although the error bars are increas¬ 
ing, the values of the densities do not change with the 
strip size. This allows us to conclude that pf*" = 0.378(1) 
and P 2 ^ = 0.2037(8) without extrapolations. As ex¬ 
pected, these values are larger than the ones found in the 
Bethe lattice solution for coordination q = A (pi Ri 0.2985 
and p 2 ~ 0.1117). 

In accordance with the solution on the Bethe lattice, 
the total density of particles for a fixed pressure - we 
may identify the reduced grand-canonical free energy fao 
in Eq. (HU) with the reduced pressure - as a function of 
the fugacity zi or of the density pi of small particles has 
a non-monotonic behavior, it displays a minimum in the 
fluid phase, as shown in Fig. [5] These curves may be 
obtained fixing the value of the dominant eigenvalue of 
the transfer matrix (this corresponds to a fixed value of 
the pressure), so that Z 2 and p may be found as func¬ 
tions of zi- Since at larger widths the finite size effects 




FIG. 5. (Color on line) (a) Transition (critical and coexis¬ 
tence) lines calculated for strips of size L in the densities’ 
space, for the fluid (full) and solid (dashed lines) phases . 
(b) Extrapolated critical (full) and coexistence (dashed) lines 
considering three point fits for the sizes {L — 2, L, L + 2) and 
the exponents indicated in the text. The square is placed at 
the localization of the tricritical point estimated by Poland 

la¬ 


in these calculations are rather small, we have done them 
for strips of width L = 18 without extrapolations. The 
locations of these minima originate a curve similar to the 
temperature of maximal density (TMD) curve delimiting 
the density anomaly in water-like fluids [l^, so we will 
call these minima as mD (minimal density) points. No¬ 
tice that in the case of water, the curves which exhibit 
a maximum are the isobars of the density as a function 
of the temperature. Here we also have isobars, but the 
density is expressed as a function of an activity, which is 
also a field-like variable in the thermodynamic sense, or 
of the conjugated density, which is a monotonic function 
of the activity. The mD curves seem to start at zi = 1/3 
as Z 2 —>■ 0 and then they are increasing functions of zi. 
Although it is very difficult to calculate the minima close 
to the tricritical point, it seems that the curves end ex¬ 
actly at this point. For small zi, the mD’s obtained for 





















FIG. 6. (Color on line) Total density of particles p = pi + 2p2, 
calculated on a strip of size L = 18, against z\ for several 
values of pressure P. From the left to the right, curves cor¬ 
respond to increased pressures. The minima in the densities 
curves define the mD line (dashed). 


different strip sizes have negligible corrections, but close 
to the tricritical point they present an appreciable L- 
dependence. Thus, we estimated the mD curve in the 
thermodynamic limit using a three-point fit. 


D. Phase diagrams 


Figure [7] summarizes our results, showing the extrap¬ 
olated critical line, coexistence line and mD curve in the 
reduced activity fraction variables Zi/{1 + Zi) in the main 
plot and activities zi, Z 2 in the inset. The phase dia¬ 
gram in the density variables is shown in Fig. [7 )d. Al¬ 
though these phase diagrams are qualitatively similar to 
the ones we obtained previously in the Bethe lattice so¬ 
lution of the model [1^, as expected the transition lines 
are shifted to larger activities (or densities) on the square 
lattice. Also, we notice that the lines of the densities of 
coexisting phases meet at an angle in the Bethe lattice 
solution, while in our results for the square lattice they 
meet tangentially. This is expected, since the classical 
value of the tricritical exponent associated to the behav¬ 
ior of the order parameter p is /32 = 1 and the estimates 
for this exponent in two dimensions, one below the upper 
tricritical dimension d = 3, is smaller than 1 231. 


IV. CONCLUSION 

In this paper we study a model of a mixed lattice gas 
with two kinds of particles. Small particles exclude only 
the site they occupy, while large particles exclude, be¬ 
sides their site, its four (square lattice) first neighbors 
also. The model with only large particles is well studied 



Z/(1+Zj) 



FIG. 7. (Golor on line) (a) Phase diagram in the variables 
Zi/(1 + Zi) (main plot) and Zi (inset), with i = 1, 2. Critical, 
coexistence and mD lines are indicated by continuous (red), 
dashed (blue) and dotted (black) lines, respectively, (b) Phase 
diagram in the densities pi and p 2 space. The densities of 
fluid (blue) and solid (green) phases at the same point of the 
coexistence line are connected by tie lines. 


in the literature and displays a continuous transi¬ 
tion between a fluid phase where sites of the lattice are 
occupied at random to a solid phase where one of the 
two sub-lattices is preferentially occupied by the parti¬ 
cles. As the small particles are introduced, a line of con¬ 
tinuous transitions is found, which ends at a tricritical 
point, beyond which the transition is discontinuous. Us¬ 
ing transfer matrix and finite-size scaling techniques, we 
estimate the thermodynamical behavior of the model ex¬ 
trapolating data of its solutions on strips of finite widths 
to the two-dimensional limit. 

The phase diagrams we found both in activity vari¬ 
ables and densities spaces are qualitatively similar to the 
ones we found before in the Bethe lattice solution [T3| . 
but some characteristics close to the tricritical point are 
distinct, since at two dimensions non-classical tricritical 
exponents are expected, leading to changes as compared 
to the mean-field behavior, as was discussed above. It is 
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interesting to compare our results with the ones obtained 
by Poland using high deirsity series expansions M- Our 
estimate for the density of small particles at the TC point 
(pf*" = 0.378(1)) is higher than the value obtained by 
Poland = 0.29 ± 0.02), but our estimate for the 

density of large particles at this point (pj*" = 0.2037(8)) 
is consistent with the one he found (pj*^ = 0.20 ± 0.01). 
We notice that we obtained estimates for the localiza¬ 
tion of the tricritical point with different methods, which 
lead to consistent results. In particular, the estimate by 
Poland is not close to our results for the critical line in 
the density variables, it is below all curves shown in Fig. 
[5l where it is represented by the black square. Also, the 
critical line in the phase diagram in the activity variables. 
Fig. [^, has a negative slope at small values of zi, which 
becomes positive after a minimum. Thus, the same re¬ 
entrant behavior found in the Bethe lattice solution of 
the model 1J| is also found in our results for the square 
lattice. This shows that, when their density is low, the 
small particles facilitate the ordering of the large ones - 
namely, they acts as an effective entropic attractive force 
among the large particles - regardless the structure of the 
underlying lattice. It may be mentioned that in a study 


of the isotropic-nematic transition for polydisperse rods, 
the critical curve also has a non-monotonic behavior in 
the space of two fugacities, although the mean length of 
the rods changes monotonically along the curve 241. 

Our estimates for the central charge are consistent with 
the Ising universality class on the whole critical line, in 
accordance with the very precise results obtained with 
similar techniques by Guo and Blbte for the model at 
the particular point of the critical line with large parti¬ 
cles only 0 and with results of a similar model which 
may be mapped, in a particular case, on the Ising model 
[lll |. The estimates for the central charge at the tricrit¬ 
ical point are close to the value of the tricritical Ising 
(BEG) universality class. 
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